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Thestandard model of particle physics describes the vast majority of experiments 
andobservations involving elementary particles. Any deviation from its predictions 
would beasign of new, fundamental physics. One long-standing discrepancy 
concerns the anomalous magnetic moment ofthe muon, a measure of the magnetic 
field surrounding that particle. Standard-model predictions! exhibit disagreement 
with measurements? that is tightly scattered around 3.7 standard deviations. Today, 
theoretical and measurement errors are comparable; however, ongoing and planned 
experiments aim to reduce the measurement error by a factor of four. Theoretically, 
the dominant source of error is the leading-order hadronic vacuum polarization 
(LO-HVP) contribution. For the upcoming measurements, it is essential to evaluate 
the prediction for this contribution with independent methods and to reduce its 
uncertainties. The most precise, model-independent determinations so far rely on 
dispersive techniques, combined with measurements of the cross-section of 
electron-positron annihilation into hadrons? 5. To eliminate our reliance on these 
experiments, here we use ab initio quantum chromodynamics (QCD) and quantum 
electrodynamics simulations to compute the LO-HVP contribution. We reach 
sufficient precision to discriminate between the measurement of the anomalous 
magnetic moment of the muon and the predictions of dispersive methods. Our result 
favours the experimentally measured value over those obtained using the dispersion 
relation. Moreover, the methods used and developed in this work will enable further 


increased precision as more powerful computers become available. 


The muon is an ephemeral sibling of the electron. It is 207 times more 
massive, but has the same electric charge and spin. Similarly to the 
electron, it behaves like a tiny magnet, characterized by a magnetic 
moment. This quantity is proportional to the spin and charge of the 
muon and inversely proportional to twice its mass. Dirac's relativistic 
quantum mechanics predicts that the constant of proportionality, g,, 
should be equal to 2. However, ina relativistic quantum field theory 
suchas the standard model, this prediction receives small corrections 
due to quantum vacuum fluctuations. These corrections are called the 
anomalous magnetic moment and are quantified by (g, - 2)/2. They 
were measured to a precision of 0.54 ppmatthe Brookhaven National 
Laboratory in the early 2000s’, and have been calculated with a com- 
parable precision (see ref." for a recent review). 

Atthislevel of precision, all ofthe interactions ofthe standard model 
contribute. The leading contributions are electromagnetic and 
described by quantum electrodynamics (QED), but the onethat dom- 
inates the theoretical error is induced by the strong interaction and 


requires solving the highly nonlinear equations of QCD at low energies. 
This contribution is determined by the LO-HVP, which describes how 
the propagation of a virtual photon is modified by the presence of 
quark and gluon fluctuations in the vacuum. Here we compute this 
LO-HVP contribution to (g, - 2)/2, denoted by aj", using ab initio 
simulations in QCD and QED. 

QCDisageneralized version of QED. The Euclidean Lagrangian forthis 
theoryis£ = [V(4e?)]F,F,,* [Vg ?Itr(G,, G.,)* Ly Gly, (0,* iqpA, + iB,)+ m, luy, 
where y, are the Dirac matrices, f runs over the ‘flavours’ of quarks, m; 
are the masses and the q;are the charges of quarks in units of the elec- 
tron charge, e. Moreover, F, = 0,A, - 0,A, and G,,, = 0,B, - 0,B,,+ [B,,, By] 
and gis the QCD coupling constant. In electrodynamics the gauge 
potential A, is a real-valued field, whereas in QCD B, is a3 x 3 Hermitian 
matrix field. The different flavours of quarks are represented by 
independent fermionic fields, y; These fields have an additional 
‘colour’ index in QCD, which runs from 1 to 3. In this work, we include 
both QED and QCD, as well as four non-degenerate quark flavours 
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Fig. 1| Contributions to a,, including examples of the corresponding 
Feynman diagrams. Solid lines are quarks and curly lines are photons. Gluons 
are not shown explicitly, and internal quark loops are shown only if they are 
attached to photons. Dots represent coordinates in position space, boxes 
denote the mass insertion relevant for strong-isospin symmetry breaking. 
The numbers give our result for each contribution; they correspond to our 


(up, down, strange and charm), in a lattice formulation that takes 
into account all dynamical effects. We also consider the tiny contribu- 
tions of the bottom and top quarks, as discussed in Supplementary 
Information. 

We compute a in the so-called time-momentum representa- 
tion?, which relies on the following two-point function with zero 
three-momentum in Euclidean time t: 


LO-HVP ; 


60-2, Y faxy, 0), (1) 


3e? H=1,2,3 


where J, i is the quark electromagnetic current, with — L = Fay, u- 
34), d- 38Y, s+ Zey c.u,d,sand care the up, down, strange and. charm 
quark fields, respectively, and the angle brackets stand for the 
QCD + QED expectation value to order e). Itis convenient to decompose 
G(t) into light, strange, charm and disconnected components, which 
have very different statistical and systematic uncertainties. Integrating 
the one-photon-irreducible part of the two-point function (equa- 
tion (1)), Giy, yields the LO-HVP contribution to the magnetic moment 
ofthe muon? *': 


gio P = gy? f , KOGO, (2) 


with weight function 


oo 2 2 
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‘reference’ system size defined byL,.;= 6.272 fm spatial and 7,47 9.408 fm 
temporal lattice extents. We also explicitly compute the finite-size corrections 
that must be added to these results, which are given separately in the lower 
right panel. The first error is the statistical and the second is the systematic 
uncertainty, except for the contributions for which only asingle, total error is 
given. Central values are medians; errors ares.e.m. 


and where w(r) = [r+ 2-J/r(r+4) /.[r(r- 4) , wis the fine-structure 
constant in the Thomson limit and m, is the muon mass. Because we 
consider only the LO-HVP contribution, for brevity we drop the super- 
script and multiply the result by 10", that is, a, stands for a!" V? x 1970 
in the following. 

The subpercent precision that we are aiming for represents a huge 
challenge for lattice QCD. To reach that goal, we must address four 
critical issues: scale determination; noise reduction; QED and strong- 
isospinsymmetry breaking; and infinite-volume and continuum extrap- 
olations. We discuss these one by one. 

The first issue is scale determination. The quantity a, depends 
on the muon mass. When computing equation (2) on the lattice, m, 
must be converted into lattice units, am,, where a is the lattice spac- 
ing. A relative error of the lattice spacing propagates into about a 
twice-as-large relative error on a,, so that a must be determined with a 
precision of few parts per thousand. We use the mass ofthe Q baryon, 
M,=1,672.45(29) MeV, from ref. ' to set the lattice spacing, where the 
uncertainty in the parentheses denotes one standard deviation. We 
also use a scale based on the gradient flow from ref. ?, denoted as Wo, 
to define an isospin decomposition of our observables. Although wo 
can be determined with sub-per-thousand precision on the lattice, it 
is inaccessible experimentally. In this work we determine the physical 
value of w, by including QED and strong-isospin symmetry-breaking 
effects: Wo = 0.17236(29)...(63),, (70). f m, where the first error is 
statistical, the second is systematic and the third is the total error. 
In total, we reach a relative accuracy of 4%o, which is better than the 
error of the previous best determination”, the value of which agrees 
with ours. There the pion decay constant was used as experimental 
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Fig. 2| Continuum extrapolation of the light connected component ofa,, 
alien, Before extrapolation we apply ataste-improvement procedure onthe 
correlator, starting at some distance t,.,. (See Supplementary Information for 
details on the improvement ‘SRHO’.) Datasets are shown for two choices of t,.,, 
0.4 fm (red) and 1.3 fm (blue). The corresponding lines show fits using linear 
and quadratic terms of à? with varying number of lattice spacings in the fit. Our 
final analysis involves about 500,000 different continuum extrapolations, 
shown inthe histogram on the left. The purple line in the left panel shows the 
central value of the final result. To estimate the error related to the 
taste-improvement procedure, we use next-to-next-to-leading-order 
staggered chiral perturbation theory (NNLO) in the long-distance part of the 
correlator (t>1.3 fm). The corresponding data are shown with grey points, 
together witha histogram, from which the systematic error related to the taste 
improvement is obtained. The total error of the final result is given by the grey 
band inthe left panel. Central values are medians; errors are s.e.m. The results 
are obtained on lattices of sizes L = 6 fm. 


input, and the isospin-symmetry-breaking effects were included only 
as an estimate. 

The second issue is noise reduction. Our result for a, is obtained 
as an integral over the conserved current-current correlation func- 
tion, from zero to infinite time separation, as shown in equation (2). 
For large separations the correlator is noisy, and this noise manifests 
itself asa statistical error ina,. To reach the desired accuracy on a, one 
needs high precision at every step. Over 20,000 configurations were 
accumulated for our 27 ensembles on L = 6 fm lattices (L is the spatial 
extent of the lattice). In addition, we include a lattice with Z ~ 11 fm. 
The most important improvement over our earlier a, determination 
in ref.“ is the extensive use of analysis techniques that are based on the 
lowest eigenmodes of the Dirac operator; see, for example, refs. 515, 
Anaccuracy gain of aboutanorder of magnitude canbe reached using 
this technique for a, (refs. ???). 

The third issue is isospin-symmetry breaking. The precision needed 
cannot be reached with pure, isospin-symmetric QCD. Thus, we 
include QED effects and allow the up and down quarks to have differ- 
ent masses. These effects are included both in the scale determination 
and in the current-current correlators. We note that the separation 
of isospin-symmetric and isospin-symmetry-breaking contributions 
requires a convention, which we discuss in detail in Supplementary 
Information. Strong-isospin breaking is implemented by taking deriva- 
tives of QCD + QED expectation values with respect to up/down quark 
masses and computing the resulting observables on isospin-symmetric 
configurations”. We note that the first derivative of the fermionic 
determinant vanishes. We also implement derivatives with respect 
to the electric charge”. It is useful to distinguish between the electric 
charge in the fermionic determinant (e, or sea electric charge) and in 
the observables (e, or valence electric charge). The complete list of 
graphs that should be evaluated are shown in Fig. 1 with our numerical 
results for them. 
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Fig. 3 | Comparison of recent results for the LO-HPV contribution tothe 
anomalous magnetic moment of the muon. See ref.’ for arecent review. 
Green squares are lattice results: this result (filled symbol) and those of 
Gérardin et al.”, Davies et al.?, Giusti et al.**, Blum et al.” and our earlier work, 
Borsanyi et al.*. Central values are medians; error bars are s.e.m. Compared to 
Borsanyi et al.*, this work has increased the accuracy of the scale setting from 
the per cent to the per thousand level; has decreased the statistical error from 
7.5 to 2.3; has computed all isospin-symmetry-breaking contributions, as 
opposed to estimating it, with the corresponding error being 1.4, down from 
5.1; has made a dedicated finite-size study to decrease the finite-size error from 
13.5to2.5; has decreased the continuum extrapolation error from 8.0 to 4.1 by 
obtaining much more statistics on our finest lattice and applying taste 
improvement. Redcircles were obtained using the R-ratio method by Davier 
etal.?, Keshavarzi et al.*, and Colangelo et al.° and Hoferichter et al.°; these 
results use the same experimental data as input. The blue shaded region is the 
value that a°"? should have to explain the experimental measurement of 


H 
(g,- 2), assuming no new physics. 


The final observable is given as a Taylor expansion around the 
isospin-symmetric, physical-mass point with zero sea and valence 
charges. Instead of the quark masses, we use the pseudoscalar meson 
masses of pions and kaons, which can be determined with high preci- 
sion. Using the expansion coefficients, we extrapolate in the charges, 
inthe strong-isospin symmetry-breaking parameter and in the lattice 
spacing, and interpolate in the quark masses to the physical point. Thus, 
we obtain a, and its statistical and systematic uncertainties. 

The fourth issue is the extrapolation to the infinite-volume and con- 
tinuum limit. The standard wisdom for lattice calculations is that M,L > 4 
should be taken, where M, is the mass of the pion. Unfortunately, this 
is not satisfactory in the present case: a, is far more sensitive to L than 
other quantities, such as hadron masses, and large volumes are needed 
to reach per-thousand accuracy. For less volume-sensitive quantities, 
we use well established results to determine the finite-volume correc- 
tions on the pion decay constant? and on charged hadron masses” *®. 
Leading-order chiral perturbation theory” and two-loop, partially 
quenched chiral perturbation theory???? for a, help to describe 
finite-size corrections, but the non-perturbative, leading-order, large-L 
expansion of ref. ? indicates that those approaches still lead to sys- 
tematic effects that are larger than the accuracy that we are aiming 
for. In addition to the infinite-volume extrapolation, the continuum 
extrapolation is also difficult. This is connected to the taste-symmetry 
breaking of staggered fermions, which we use in this work. 

We correct for finite-volume effects on a, by computing them directly 
by performing lattice simulations on L « 11 fm lattices, with highly 
suppressed taste violations and with physical, taste-averaged pion 
masses. These corrections are cross-checked against three models 
that describe the relevant long-distance physics, in turn validating 
the use of these models for the residual, sub-per-thousand extrapola- 
tion to infinite volume. These models include: (i) the full two-loop, 
finite-volume, chiral perturbation theory corrections for a,; (ii) the 
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Fig. 4 | Continuum extrapolation ofthe isospin-symmetric, light, 
connected component ofthe window observable a, win, (aie iso The data 
points are extrapolated to the infinite-volume limit. Central values are 
medians; error bars ares.e.m. Two different ways to perform the continuum 
extrapolations are shown: one without improvement, and another with 
corrections from a model involving the p meson (SRHO). In both cases the lines 
show linear, quadratic and cubic fits in à? with varying number of lattice 
spacings in the fit. The continuum-extrapolated result is shown with the results 
from Blum etal.” and Aubin et al.”. Also plotted is our R-ratio-based 
determination, obtained using the experimental data compiled by the authors 
of ref.* and our lattice results for the non-light-connected contributions. This 
plotis convenient for comparing different lattice results. Regarding the total 
Q, wi; for which we must also include the contributions of flavours other than 
light and isospin-symmetry-breaking effects, we obtain 236.7(1.4),,, onthe 
lattice and 229.7(1.3),: from the R-ratio; the latter is 3.70 or 3.1% smaller than the 
lattice result. 


Meyer-Lellouch-Lüscher-Gounaris-Sakurai technique described in 
Supplementary Information; and (iii). the p-7r-y model ofJegerlehner 
and Szafron”, already used in a lattice context in ref. ?'. Moreover, to 
reduce discretization errors in the light-quark contributions to a,, 
before extrapolating those contributions to the continuum, we apply 
ataste-improvement procedure that reduces lattice artefacts due to 
taste-symmetry breaking. The procedure is built upon the three models 
of m-p physics mentioned above. We provide evidence that validates 
this procedure in Supplementary Information. 

Combining all of these ingredients, we obtain as a final result 
Ay = 707.5(2.3) sat(5-O)syst(5-5) ror The statistical error comes mainly 
fromthe noisy, large-distance region ofthe current-current correla- 
tor. Thesystematic erroris dominated by the continuum extrapola- 
tion and the finite-size effect computation. The total error is obtained 
by adding the first two in quadrature. In total, we reach a relative 
accuracy of 0.8%. In Fig. 2 we show the continuum extrapolation of 
the light, connected component of a,, which gives the dominant 
contribution to q,. 

Figure3 compares our result with previous lattice computations and 
also with results from the R-ratio method, which have recently been 
reviewed in ref. ". In principle, one can reduce the uncertainty of our 
result by combining our lattice correlator, G(t), with the one obtained 
fromthe R-ratio method, in regions of Euclidean time in which the lat- 
ter is more precise”. We do not do so here because there is a tension 
between our result and those obtained by the R-ratio method, as canbe 
seeninFig.3.Forthetotal LO-HVP contribution to a, our result is 2.06, 
2.50,2.40and 2.20 larger than the R-ratio results of a, = 694.0(4.0) (ref.?), 
à, = 692.78(2.42) (ref. *), a, = 692.3(3.3) (refs. °°) and the combined 
result a, = 693.1(4.0) of ref. ’, respectively. It is worth noting that the 
R-ratio determinations are based on the same experimental datasets 
andare therefore strongly correlated, although these datasets were 
obtained in several different and independent experiments that we have 


54 | Nature | Vol 593 | 6 May 2021 


noreasonto believe are collectively biased. Clearly, these comparisons 
need further investigation, although it should also be kept in mind 
that the tensions observed here are smaller, for instance, than what 
is usually considered experimental evidence for a new phenomenon 
(30) and much smaller than what is needed to claim an experimental 
discovery (5o). 

As afirst step in that direction, it is instructive to consider a mod- 
ified observable, where the correlator G(f) is restricted to a finite 
interval by asmooth window function”. This observable, which we 
denote as A,,win, is obtained much more readily than a, on the lattice. 
Its shorter-distance nature makes it far less susceptible to statistical 
noise and to finite-volume effects. Moreover, in the case of staggered 
fermions, it has reduced discretization artefacts. This is shown in 
Fig. 4, where the light, connected component of a, is plotted as 
a function of a°. Because the determination of this quantity does 
not require overcoming many of the challenges described above, 
other lattice groups have obtained it with errors comparable to 
ours!?”°, This allows a sharper benchmarking of our calculation of 
this challenging, light-quark contribution that dominates a,. 
Our a! differs by 0.20 and 2.20 from the lattice results of ref. ?? 
and ref. ?, respectively. Moreover, a, yi, can be computed using the 
R-ratio approach, and we do so using the dataset provided by the 
authors of ref. *. However, here we find a 3.70 tension with our lattice 
result. 

To conclude, when combined with the other standard-model con- 
tributions (see, for example, refs. **), our result for the leading-order 
hadronic contribution to the anomalous magnetic moment of the 
muon, d," = 707.5(5.5) ror x 10? , weakens the long-standing dis- 
crepancy between experiment and theory. However, as discussed above 
and can be seen in Fig. 2, our lattice result shows some tension with the 
R-ratio determinations of refs. *°. Obviously, our findings should be 
confirmed-or refuted—by other studies using different discretizations 
of QCD. Those investigations are underway. 
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Methods 


Finite-size effects 

Finite-size effects on a, were the largest source of uncertainty in our 
previous work'*. We compute these effects in a systematic way, which 
includes dedicated lattice simulations, chiral perturbation theory 
(XPT) and phenomenological models. The goal is to provide a single 
number that is to be added to the continuum-extrapolated lattice 
result obtained in a reference box, which is defined by a spatial extent 
of L,.¢= 6.272 fm and a temporal extent of T, = (2/2)L,,.. Here we sum- 
marize our findings on the finite-size effect of the isospin-symmetric 
part. More details anda discussion oftheisospin-symmetry-breaking 
part can be found in Supplementary Information. 

We perform dedicated lattice simulations with two different lat- 
tice geometries: oneisa56 x 84 lattice with the reference box size and 
the other is a large 96 x 96 lattice with box size L = Ly, = 10.752 fm and 
T= Ty 7 Ly, Because taste violations distort the finite-size effects, we 
designed a new action with highly suppressed taste breaking, which 
we call 4HEX. Our strategy is then to compute the finite-size correction 
as the following sum: 


a,(~, oo) — Ay (Lret, Tref) 
7[a,(Loig, Toig) 7 U(Lret Tref) lanex (4) 
+[a,,(%, *) - a,(Lyis, Toig) Ixpr- 


The first difference on the right-hand side is taken from the dedi- 
cated 4HEX simulations. The second difference is expected to be 
much smaller than the first, and is taken from a non-lattice approach: 
two-loop XPT. 

We consider four non-lattice approaches to compute both differ- 
ences on the right-hand side of equation (4). In the case of the first 
difference, the results obtained are compared to our 4HEX simula- 
tions. The first approach is XPT to next-to-leading order (NLO) and 
next-to-next-to-leading order (NNLO), the second is the Meyer-Lel- 
louch-Lüscher-Gounaris-Sakurai (MLLGS) model, the third approach 
is that of Hansen and Patella (HP)? and the fourth is the p-7-y (RHO) 
model of ref. ”. 

We compute the first difference in equation (4) using dedicated 
simulations with the 4HEX action. We use the harmonic mean square 
(HMS) to set the physical point: 


z 1 cE 
Mums = 16 2 Mi), 
a 


defined as an average over the masses of the 16 pion tastes, M, ,. We set 
M,,ums to the physical value of the pion mass, which requires lowering 
the Goldstone pion mass to 110 MeV. This way of fixing the physical 
point results in much smaller lattice artefacts than the usual setting 
with the Goldstone pion, at least for an observable suchas the finite-size 
effect. 

To generate the 4HEX dataset, we performed simulations with two 
different Goldstone pion masses: M,,=104 MeV and 121 MeV. To set the 
physical point as described above, we performed an interpolation from 
these two pion masses to M, - 110 MeV. 

To compute ae from the current propagator in our 4HEX simula- 
tions, we use the upper- and lower-bound technique described in Sup- 
plementary Information. Results for the M,=121 MeV simulation point 
are plotted in Extended Data Fig. 1. The bounds meet at around 4.2 fm 
and 4.7 fm on the small and large volumes, respectively. At these dis- 
tances, we take the average of the two bounds as an estimate for alent, 
The results are given in Extended Data Table 1. 

We have only one lattice spacing with the 4HEX action, so the 
finite-size effects cannot be extrapolated to the continuum limit. We 
estimate the cutoff effect of the result by comparing the total a, with the 


4HEX action atthis single lattice spacing to the continuum-extrapolated 
4stout action-based lattice result, both in the L,e volume. The 4HEX 
result is about 7% larger than the continuum value. Therefore, we reduce 
the measured finite-size effect by 7?6, and assign a 7?6 uncertainty to 
this correction step. For the difference we get: 


Gy Lyig, Toig) a A (Ler, Teen) = 18.1(2.0) sear (1-4) cont: (5) 


The result is obtained from the ane numbers of Extended Data 
Table 1, including a (9/10) charge factor. The first error is statistical and 
the second is an estimate of the cutoff effect. 

Extended Data Table 2 collects the finite-size effect computed in 
various non-lattice approaches. The different models give finite-size 
effects of similar size, which agree well with the lattice determination 
of equation (5). Only the NLO result differs by about 3o. The fact that 
NLO XPT underestimates the finite-size effect was shown in ref. 5, at 
anon-physical pion mass. Using the physical pion mass, a dedicated 
finite-volume study was presented in ref. °°; it reaches the same conclu- 
sion as we do, albeit with larger errors. We also see that, according to 
the models, the finite-T effect is much smaller than the finite-L effect. 

The good agreement for the finite-size effect of the reference box 
between the models and the lattice gives us confidence that the models 
can be used toreliably compute the very small, residual, finite-size 
effect of the large box. The corresponding model estimates can be 
found in Extended Data Table 1. For an infinite-time extent, the NNLO 
XPT, HP and RHO approaches agree nicely. As a final value for the 
large-box, finite-size effect we take the NNLO XPT result including 
finite-7 effects: 


ae, œ) — a, (Lpig, Toig) = 0.6(0.3)pig, 


where the uncertainty is an estimate of higher-order effects, given here 
by the difference of the NNLO and NLO values. 

For our final result for the finite-size effect of the reference 
box, we also include the contribution of isoscalar channel and 
isospin-symmetry-breaking effects, giving: 


ay, eo) - GLrer, Tref) 
= 18.7(2.0) rat (1-4) cont(0-3)pig(O-6)-0(0-1) gea(2-5) ot 


The first error is the statistical uncertainty of our 4HEX computa- 
tion, the second is an estimate of the 4HEX cutoff effects, the third is 
the uncertainty of the residual finite-size effect of the ‘big’ lattice, the 
fourth is an XPT estimate of the isoscalar (isospin /= 0) finite-size effect 
and the fifth is an estimate of the isospin-symmetry-breaking effects. 
The last, total error is the sum of the first five, added in quadrature. 
The vast majority of the finite-size effects is obtained using the 4HEX 
lattice computation; for the rest we apply analytical methods that are 
validated by the lattice computation: for the main contribution they 
give values that are consistent with the lattice result. 


Taste improvement 
As is well known, some of the most important cutoff effects of stag- 
gered fermions are taste violations. At long distances, these violations 
distort the pion spectrum. Because a, is predominantly a long-distance 
observable, dominated by a two-pion contribution, including the p 
resonance, we expect these effects to be largest in the light-quark terms. 
Weinvestigate various physically motivated models for reducing 
long-distance taste violations in our lattice results. We consider three 
techniques: NNLO XPT, the MLLGS model and the RHO model. For 
the definition of these models, see Supplementary Information. We 
investigate and discuss the suitability of their staggered versions for 
reducing the taste violations present in our lattice data. We call the 
resulting corrections taste improvements, because they improve the 
continuum extrapolation of our lattice data without, in principle, 


modifying the continuum-limit value. Indeed, these corrections vanish 
inthat limit, as taste-symmetry-breaking effects should. These improve- 
ments are applied on light-quark observables at the isospin-symmetric 
point, the taste violations of which have the largest impact on our final 
uncertainties. 

The NNLO XPT, MLLGS and RHO models describe the long-distance 
physics associated with finite-volume effects, as measured in our 
simulations. One can also define corresponding models that describe 
thetaste violations, denoted as NNLO SXPT, SMLLGS and SRHO, respec- 
tively. We find that these describe the physics associated with taste 
violations, at least at larger distances. This is illustrated in Extended 
Data Fig. 2, where cutoff effects in the integrand of a2" are plotted as 
afunction of Euclidean time. More specifically, we define the physical 
observable obtained by convoluting the integrand of ash witha 
smooth window function W(t; t) of a width of 0.5 fm and starting at a 
time of tj. Then we consider the difference in the value of this observ- 
able, obtained on a fine and a coarse lattice at a sequence of t; values 
separated by O.1fm. These are compared to the NLO SXPT, NNLO SXPT, 
SRHO and SMLLGS predictions for this quantity, evaluated at the exact 
parameters ofthe ensembles. 

The SMLLGS, SRHO and NNLO SXPT taste improvements describe 
the numerical data very nicely for t, z 2.0 fm, fairly well for t, z 1.0 fm 
and all the way down to t, = 0.4 fm in the case of SRHO. All three slightly 
overestimate the observed cutoff effects, with the p-meson-based 
approach performing best, whereas NNLO displays a large deviation 
from the lattice results in the t, x 0.8 fm region. The lattice results have 
a maximum at ¢,=1.4 fm, as does the SRHO improvement, reinforcing 
our confidence that this model captures the relevant physics. 

These findings lead us to apply the following taste corrections to 
our simulation results for alg, T, a), obtained on an Z? x T lattice 
with lattice spacing a, before performing continuum extrapolations: 


ale" (I, T, a) > a£" (L, T, a) 


+ OL aio, Los Tre) - a8 (LT, a) 

with tep = 0.4, 0.7, 1.0 and 1.3 f m, and where the factor (10/9) is related 
to the quark charges. We note that by using L,e and 7,,,in the above 
equation, we apply a very small volume correction to interpolate all of 
our simulation results to the same reference four-volume, so that they 
can be extrapolated to the continuum limit together. 

The taste-improved data are then extrapolated to the contin- 
uum using our standard fitting procedure, in the course of which 
isospin-symmetry-breaking effects are also included. To estimate the 
systematic error we use a histogram technique. The central values and 
the detailed error budget of this analysis can be found in Supplemen- 
tary Information. 

The procedure described above does not yet take into account the 
systematic uncertainty associated with our choice of SRHO for taste 
improvement for t > 1.3 fm. Given that applying no taste improve- 
ment in that region is not an option, because of the nonlinearities 
introduced by two-pion taste violations, we turn to NNLO SXPT, only 
as ameans of estimating the uncertainty associated with this choice. 
Thus, we define this systematic uncertainty as ERR = (SRHO - NNLO 
SXPT) for t» 1.3 f m. Then, we perform the same histogram analysis 
but with SRHO, SRHO - ERR and SRHO + ERR improvements. From 
this histogram we extract the contribution that comes from the vari- 
ation inthe improvement model from SRHO - ERR to SRHO + ERR. We 
assign this full spread to the systematic uncertainty associated with 
the taste-improvement procedure. We add this error in quadrature to 
the error given by the histogram technique discussed in the previous 
paragraph. 

The procedure is illustrated in Extended Data Fig. 3, which shows 
the datasets for a'£"* without and with taste improvements, as func- 
tions of a. (See also Fig. 2, which zooms in on the taste-improved, 


continuum extrapolations.) The SRHO improvements with £,., = 0.4 fm 
are shownas red points, whereas blue points correspond to ¢,.,=1.3 fm. 
These plots include isospin-symmetry-breaking contributions. An 
example of our lattice results with SRHO improvement between 
t=0.4 fm and t=1.3 fmand NNLO SXPT improvement are shown as grey 
points in Extended Data Fig. 3. 

Animportant check of our taste-improvement procedure is provided 
by the study of the isoscalar or /= O contribution to a,, as suggested 
by arguments made in ref. ?. Here we work with the isospin-symmetric 
datasets. Thenthe/- O contribution is defined as: 


e i 1 
q^? = —_qiight ab ads + ue 


n Io% i (6) 


where the ellipsis stands for the quark-connected contributions of the 
more massive ss, c,... quarks. This quantity receives no two-pion contri- 
butions: it starts with three pions, thetaste-symmetry-breaking effects 
of which should be very small. Thus, if our understanding of discretiza- 
tion errors in aj" and a?** is correct, the taste-symmetry-breaking 
corrections observed in the light and disconnected quantities must 
belargely absent from d. Asaconsequence, we expect the continuum 
extrapolation ofa; ? to be much milder. That is exactly what is Shown 


in Extended Data Fig. 4 and explained in its caption. 


Results for wy, M,,and AM 

Here we describe the details of the analyses that are used to obtain the 
physical values of Wo, M,, and AM? = M2, - M2, from the experimental 
values of hadron masses, including the mass of the Q baryon. M,,, M 
and M,, denote the masses of mesons built from an up and an anti-up, 
a down and an anti-down, and a strange and an anti-strange quark, 
respectively, without quark-disconnected contributions. The results 
for these quantities are used to define the isopin-symmetric point in 
this work, as described in Supplementary Information. 

In the case of Wo, the observable that we fit is wWoMo. To account for 
the systematic error due to the different continuum extrapolations, 
we apply both linear and quadratic functions in the isospin-symmetric 
component. We also remove zero/one/two/three of the coarsest lat- 
tice spacings in the linear fit and zero/one/two lattice spacings in the 
quadratic fit. For the tiny valence QED component only linear fits are 
applied, with zero/one/two points removed; for the even smaller sea 
QED contributions, we have either a constant or a linear fit with all 
lattice spacings. 

The systematic error of the hadron mass fits is taken into account 
by 24 different combinations of the fit ranges: three for the M; mass, 
two for the pseudoscalars, two for the isospin breaking of the M,, and 
two for the isospin breaking of the pseudoscalars. To account for the 
experimental error on Mo, we carry out the analysis with two different 
experimental values: one that corresponds to the central value plus the 
experimental error, and another with this error subtracted. 

Altogether, these yield a total of 129,024 fits. When the different 
analyses are combined into a histogram to determine the systematic 
error, the results from different fitting functions or lattice spacing cuts 
are weighted with the Akaike information criterion and the rest with 
flat weighting. We obtain: 


Wy = 0.17236(29) seat (63)syse(7OVeor fm, (7) 


where the first error is statistical, the second is systematic and the third 
is the total error; we reacha relative precision of 0.4%. The splitting up 
of the error into different sources is explained in Supplementary Infor- 
mation. In Extended Data Fig. 5 we show the various isospin components 
of w,M, versus the lattice spacing squared, together with the different 
continuum extrapolations. Our result (equation (7)) is in good 
agreement with earlier four-flavour determinations: Wo = 0.1715(9) fm 
from ref. ? and w= 0.171433 fm from ref. ””. In those studies, the 
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isospin-symmetry-breaking effects were only estimated, whereas in 
our case they are fully accounted for. 

The same procedure is used for M,,. We work with (M,,/M,)” instead 
of M,JM,, because the fit qualities are much better in the first case. The 
129,024 different fits give: 


Mss = 689.89(28) sear(40)syse(49)ror MeV. (8) 


Finally, we carry out the analysis for AM?/M2, with AM? = M2, - M2,, 
whichis a measure of the strong-isospin breaking. Altogether we have 
3,328 fits, which give the following central value with statistical, sys- 
tematic and total errors: 


AM? = 13,170(320) tat(270)ys¢(420) ror MeV?. (9) 


Data availability 


The datasets used for the figures and tables are available from the cor- 
responding author on request. 


Code availability 


ACPU code for configuration production and measurements can be 
obtained from the corresponding author upon request. The Wilson flow 
evolution code, which was used to determine wọ, can be downloaded 
from https://arxiv.org/abs/1203.4469. 
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Extended Data Fig. 1| Upper and lower bounds onthe light 
isospin-symmetric component ofa,, fais". The bounds are computed 
using the lattice correlator below atime separation of t, and an analytical 
formula describing the large-time behaviour above ¢,.The results shown are 
obtained with the 4HEX action on two different lattice sizes, 56 x 84 and 


96 x 96, bothat a=0.112 fm lattice spacing and M, - 121 MeV Goldstone pion 
mass. We also carried out another simulation with M, - 104 MeV mass. From 
these two, we interpolate to M, - 110 MeV. This value ensures that a particular 


average of pion tastes is fixed to the physical value of the pion mass (see text). 
Error bars are statistical errors (s.e.m.). 
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ExtendedDataFig.2|Isospin-symmetric component of alist, computed fromthe simulation, and errors are statistical (s.e.m.). The coloured curves are 
witha sliding window. The window starts at t, and ends 0.5 fm later. The plot the predictions of NLO,NNLO SXPT, and the SRHO and SMLLGS models. They 
shows the difference between a fine and a coarse lattice with spacing are computed at the parameters (pion mass, taste violation, volume) of the 


a=0.064 fm and a=0.119 fm. The black squares with error bars are obtained simulations. 
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improvement. The blue squares repesent data that have undergone no taste SRHO improvement, and the corresponding central value and error is the 
improvement for ¢<1.3 fm and SRHO improvement above. The blue curves purple band. The darker grey circles correspond to results corrected with 
correspond to example continuum extrapolations of improved data to SRHO in the range 0.4-1.3 fm and with NNLO SXPT for larger t. These latter fits 
polynomials in @?, up to and including a*. We note that extrapolations in serve to estimate the systematic uncertainty of the SRHO improvement. The 
@a,,(1/a)?, with a,(1/a) the strong coupling at the lattice scale, are also grey band includes this uncertainty, and the corresponding histogram is shown 


considered in our final result. The red circles and curves are the same as the with grey. Errors ares.e.m. 
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Extended DataFig.4 | Comparison ofthe continuum extrapolation of 

a1, ^? to those of ali" and a‘. Top, grey points correspond to our 
uncorrected results for wae. The red symbols show the same results with our 
standard SRHO taste improvement. They havea much milder continuum limit 
thatexhibits none ofthe nonlinear behaviour ofthe grey points. The red curves 
show typical examples of illustrative continuum extrapolations of those 
points. Bottom, grey and red points and curves are the same quantities, but for 
d. Combining the results from the two individual continuum extrapolations 
of ial and a‘'*, according to equation (6), gives the result with statistical 
errors illustrated by the red band, with combined statistical and systematic 


errors indicated by the broader pink band. The blue points correspond to our 


results for aj, Orient for each of our simulations, and are obtained by combining 
the two sets of grey points, according to equation (6). As these blue points 
show, the resulting continuum-limit behaviour of aj" is much milder than that 
of either the uncorrected att or ai, and shows none of their curvature. This 
behaviour resembles much more that of the taste-improved red points. 
Moreover, all of the blue points, including typical continuum extrapolations 
drawn as blue lines, lie within the bands. This suggests that our taste 
improvements neither bias the central values of our continuum-extrapolated 
alg and ade, nor dothey lead to an underestimate of uncertainties. Errorsare 
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Extended Data Fig. 5| Continuum extrapolations ofthe contributions to 
WM). Fromtop to bottom: isospin-symmetric, electromagnetic valence- 
valence, sea-valence and sea-sea components. The results are multiplied by 
10ł/Mġ and the electric derivatives are multiplied by e2, where the asterisk 
denotes physical value. Error bars show statistical errors (s.e.m.). Dashed lines 
arecontinuum extrapolations, showing illustrative examples from our several 
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thousandfits. Only thelattice spacing dependenceis shown:the data points 
are moved to the physical light- and strange-quark mass point. This adjustment 
varies from fit to fit, and the red data points are obtained in an a?-linear fit to all 
ensembles. Ifina fit the adjusted points differ considerably fromthe red 
points, we showthem with grey colour. The final result is obtained froma 
weighted histogram of the several thousand fits. 
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Extended Data Table 1| Isospin-symmetric component of an" 


M, in 4HEX —^ | 104 MeV | 121 MeV | 110 MeV 
aisht(56 x 84) | 685.9(2.7) | 668.3(2.0) | 679.5(1.9) 
aj£" (96 x 96) | 710.7(1.9) | 684.3(1.7) | 701.1(1.3) 


ations 4HEX action on two different volumes and two different Goldstone 
physical point, using the HMS averaged p ion mass prescription. 


pion masses. In the last column we give the interpolated value at the 


Extended Data Table 2 | Finite-size effect in the reference box of the isospin-symmetric component of a,, 


NLO XPT | NNLO XPT | MLLGS | HP | RHO 
Gl Nm lcudEecg di 15.7 wm | =| = 
IS 00) — a aoo) 11.2 15.3 17.4 |163 | 148 
a, (oo oo) — as Lose, Tug) 0.3 0.6 = v 
a4 (00,00) — au Dus, 00) 1.2 1.4 — 1.4 | 1.4 


The values are obtained from various model approaches. 


